{
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.10.14",
   "mimetype": "text/x-python",
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "pygments_lexer": "ipython3",
   "nbconvert_exporter": "python",
   "file_extension": ".py"
  },
  "kaggle": {
   "accelerator": "nvidiaTeslaT4",
   "dataSources": [],
   "dockerImageVersionId": 30776,
   "isInternetEnabled": true,
   "language": "python",
   "sourceType": "notebook",
   "isGpuEnabled": true
  }
 },
 "nbformat_minor": 5,
 "nbformat": 4,
 "cells": [
  {
   "cell_type": "code",
   "source": "import torch\nfrom collections import OrderedDict\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport scipy.io\n#from pyDOE import lhs#拉丁超立方抽样\nfrom torch import nn\nimport time\nimport warnings\nwarnings.filterwarnings(\"ignore\", category=UserWarning)  # 只忽略UserWarning类型的警告",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-09-24T08:23:44.161933Z",
     "iopub.execute_input": "2024-09-24T08:23:44.162356Z",
     "iopub.status.idle": "2024-09-24T08:23:45.981338Z",
     "shell.execute_reply.started": "2024-09-24T08:23:44.162314Z",
     "shell.execute_reply": "2024-09-24T08:23:45.980437Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2025-01-06T07:50:49.832915Z",
     "start_time": "2025-01-06T07:50:49.819378Z"
    }
   },
   "id": "af4b048c0e49992d",
   "execution_count": 2,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": "if torch.cuda.is_available():\n    device = torch.device('cuda')\nelse:\n    device = torch.device('cpu')\ndef setup_seed(seed):\n    torch.manual_seed(seed)\n    torch.cuda.manual_seed_all(seed)\n    np.random.seed(seed)\n    #  random.seed(seed)\n    torch.backends.cudnn.deterministic = True\n\nsetup_seed(124)\ndtype = torch.double",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-09-24T08:23:45.986195Z",
     "iopub.execute_input": "2024-09-24T08:23:45.988477Z",
     "iopub.status.idle": "2024-09-24T08:23:46.049102Z",
     "shell.execute_reply.started": "2024-09-24T08:23:45.988430Z",
     "shell.execute_reply": "2024-09-24T08:23:46.048271Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2025-01-06T07:50:50.359320Z",
     "start_time": "2025-01-06T07:50:50.343278Z"
    }
   },
   "id": "9e2f9d35aa33fbb9",
   "execution_count": 3,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": "#求导\ndef gradients(outputs, inputs):\n    return torch.autograd.grad(outputs, inputs, grad_outputs=torch.ones_like(outputs), create_graph=True)\n#类型转换\ndef to_numpy(input):\n    if isinstance(input, torch.Tensor):\n        return input.detach().cpu().numpy()\n    elif isinstance(input, np.ndarray):\n        return input\n",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-09-24T08:23:46.052132Z",
     "iopub.execute_input": "2024-09-24T08:23:46.052502Z",
     "iopub.status.idle": "2024-09-24T08:23:46.057747Z",
     "shell.execute_reply.started": "2024-09-24T08:23:46.052456Z",
     "shell.execute_reply": "2024-09-24T08:23:46.056804Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-12-10T01:45:43.085971Z",
     "start_time": "2024-12-10T01:45:43.071903Z"
    }
   },
   "id": "e501d8818e912d51",
   "execution_count": 3,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "#构建网络layers\n",
    "def Dnnlayers(input_layers,output_layers,hidden_layers,neural):\n",
    "    layers = []\n",
    "    for i in range(hidden_layers + 2):\n",
    "        if i == 0:\n",
    "            layers.append(input_layers)\n",
    "        elif i == hidden_layers + 1:\n",
    "            layers.append(output_layers)\n",
    "        else:\n",
    "            layers.append(neural)\n",
    "    return layers\n",
    "#初值\n",
    "def IC(x, crhoL, cuL, cpL, crhoR, cuR, cpR):\n",
    "    N = x.shape[0]\n",
    "    rho_init = np.zeros((x.shape[0], 1))\n",
    "    u_init = np.zeros((x.shape[0], 1))\n",
    "    p_init = np.zeros((x.shape[0], 1))\n",
    "\n",
    "    # rho, p - initial condition\n",
    "    for i in range(N):\n",
    "        if (x[i, 1] < 0.0):\n",
    "            rho_init[i] = crhoL\n",
    "            u_init[i] = cuL\n",
    "            p_init[i] = cpL\n",
    "        else:\n",
    "            rho_init[i] = crhoR\n",
    "            u_init[i] = cuR\n",
    "            p_init[i] = cpR\n",
    "    U_ic = np.hstack([rho_init, u_init, p_init])\n",
    "\n",
    "    return U_ic\n",
    "#外推边界\n",
    "def BC(x, crhoL, cuL, cpL, crhoR, cuR, cpR):\n",
    "    N = x.shape[0]\n",
    "    rho_bc = np.zeros((x.shape[0], 1))\n",
    "    u_bc = np.zeros((x.shape[0], 1))\n",
    "    p_bc = np.zeros((x.shape[0], 1))\n",
    "\n",
    "    # rho, p - initial condition\n",
    "    for i in range(N):\n",
    "        if (x[i, 1] < 0.0):\n",
    "            rho_bc[i] = crhoL\n",
    "            u_bc[i] = cuL\n",
    "            p_bc[i] = cpL\n",
    "        else:\n",
    "            rho_bc[i] = crhoR\n",
    "            u_bc[i] = cuR\n",
    "            p_bc[i] = cpR\n",
    "    U_bc = np.hstack([rho_bc, u_bc, p_bc])\n",
    "\n",
    "    return U_bc\n"
   ],
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-09-24T08:23:46.060253Z",
     "iopub.execute_input": "2024-09-24T08:23:46.061124Z",
     "iopub.status.idle": "2024-09-24T08:23:46.074353Z",
     "shell.execute_reply.started": "2024-09-24T08:23:46.061075Z",
     "shell.execute_reply": "2024-09-24T08:23:46.073511Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-12-10T01:45:43.102065Z",
     "start_time": "2024-12-10T01:45:43.086975Z"
    }
   },
   "id": "79e6a971aae89fdf",
   "execution_count": 4,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": "class DNN(nn.Module):\n    def __init__(self, layers):\n        #继承父类\n        super(DNN, self).__init__()\n\n        #depth of network\n        self.depth = len(layers) - 1\n        #print(self.depth)\n\n        #activation of network\n        self.activation = nn.Tanh()\n\n        #create the neural network\n        layers_list = list()\n        for i in range(self.depth - 1):\n            layers_list.append(\n                ('layer_%d' % i, nn.Linear(layers[i], layers[i+1]))#create each network\n            )\n            layers_list.append(\n                ('activation_%d' % i, self.activation)  #create each activation of network\n            )\n        layers_list.append(\n            ('layer_%d' % (self.depth - 1), nn.Linear(layers[-2], layers[-1]))\n            #last network do not have activation\n        )\n        #创建一个有序字典，其中包含了从 layers_list 中获得的键值对，这在需要保持元素顺序的场景（如神经网络层的顺序）中非常有用。\n        layerDict = OrderedDict(layers_list)\n\n        #deploy layers\n        self.layers = nn.Sequential(layerDict)\n        #print(self.layers)\n\n    #forword network,output the result of network\n    def forword(self, x):\n        out = self.layers(x)\n        return out",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-09-24T08:23:46.075428Z",
     "iopub.execute_input": "2024-09-24T08:23:46.075700Z",
     "iopub.status.idle": "2024-09-24T08:23:46.086317Z",
     "shell.execute_reply.started": "2024-09-24T08:23:46.075671Z",
     "shell.execute_reply": "2024-09-24T08:23:46.085503Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-12-10T01:45:43.118014Z",
     "start_time": "2024-12-10T01:45:43.105574Z"
    }
   },
   "id": "2b5f4f99b1b8181f",
   "execution_count": 5,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": "class Pinns():\n    def __init__(self, layers):\n        self.dnn = DNN(layers).double().to(device)    \n    #forword network,output the result of network\n    def forword(self, x):\n        y = self.dnn.forword(x)\n        rho, u, p = y[:, 0:1], y[:, 1:2], y[:, 2:]\n        return rho, u, p\n    #function loss,残差点训练\n    def loss_pde(self, x_int):\n        y = self.dnn.forword(x_int)\n        rho, u, p = y[:, 0:1], y[:, 1:2], y[:, 2:]\n        #gamma = self.gamma\n\n        U2 = rho * u\n        U3 = 0.5 * rho * u ** 2 + p / (1.4 - 1.0)#0.4 #gamma - 1 = 0.4\n        # F1 = U2\n        F2 = rho * u ** 2 + p\n        F3 = u * (U3 + p)\n\n        # Gradients and partial derivatives\n        drho_g = gradients(rho, x_int)[0]\n        rho_t, rho_x = drho_g[:, :1], drho_g[:, 1:]\n\n        dU2_g = gradients(U2, x_int)[0]\n        U2_t, U2_x = dU2_g[:, :1], dU2_g[:, 1:]\n        dU3_g = gradients(U3, x_int)[0]\n        U3_t, U3_x = dU3_g[:, :1], dU3_g[:, 1:]\n        dF2_g = gradients(F2, x_int)[0]\n        F2_t, F2_x = dF2_g[:, :1], dF2_g[:, 1:]\n        dF3_g = gradients(F3, x_int)[0]\n        F3_t, F3_x = dF3_g[:, :1], dF3_g[:, 1:]\n        \n        f = (((rho_t + U2_x) ) ** 2).mean() + \\\n            (((U2_t + F2_x) ) ** 2).mean() + \\\n            (((U3_t + F3_x) ) ** 2).mean()  #\n        \n        return f\n    #initial loss function\n    def loss_ic(self, x_ic, U_ic):\n        y_ic = self.dnn.forword(x_ic)\n        rho_ic_nn, u_ic_nn, p_ic_nn = y_ic[:, 0:1], y_ic[:, 1:2], y_ic[:, 2:]\n        loss_ics = ((rho_ic_nn - U_ic[:,0:1]) ** 2).mean()+\\\n                   ((u_ic_nn - U_ic[:,1:2]) ** 2).mean()+\\\n                   ((p_ic_nn - U_ic[:,2:]) ** 2).mean()\n        return loss_ics\n    #boundary condition\n    def loss_bc(self, x_bc, U_bc):\n        y_bc = self.dnn.forword(x_bc)\n        rho_bc_nn, u_bc_nn, p_bc_nn = y_bc[:, 0:1], y_bc[:, 1:2], y_bc[:, 2:]\n        loss_bcs = ((rho_bc_nn - U_bc[:, 0:1]) ** 2).mean() + \\\n                   ((u_bc_nn - U_bc[:, 1:2]) ** 2).mean() + \\\n                   ((p_bc_nn - U_bc[:, 2:]) ** 2).mean()\n        return loss_bcs",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-09-24T08:23:46.087638Z",
     "iopub.execute_input": "2024-09-24T08:23:46.088500Z",
     "iopub.status.idle": "2024-09-24T08:23:46.103784Z",
     "shell.execute_reply.started": "2024-09-24T08:23:46.088466Z",
     "shell.execute_reply": "2024-09-24T08:23:46.102979Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-12-10T01:45:43.133193Z",
     "start_time": "2024-12-10T01:45:43.120277Z"
    }
   },
   "id": "5a835b6157600a38",
   "execution_count": 6,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "class conbine():\n",
    "    def __init__(self, layers1, layers2, layers3):\n",
    "        self.model1 = Pinns(layers1)\n",
    "        self.model2 = Pinns(layers2)\n",
    "        self.model3 = DNN(layers3).double().to(device)    \n",
    "    def loss_pde(self, x_int):\n",
    "        t = x_int[:, 0:1]\n",
    "        x = x_int[:, 1:2]\n",
    "        x_predict = self.model3.forword(t)[:,0:1]\n",
    "        \n",
    "        index1 = [i for i in range(len(x_int)) if x[i] < x_predict[i]]\n",
    "        index2 = [i for i in range(len(x_int)) if x[i] >= x_predict[i]]\n",
    "        x1_int = x_int[index1]\n",
    "        x2_int = x_int[index2]\n",
    "        #loss_pde(self, x_int):\n",
    "        loss_pde1 = self.model1.loss_pde(x1_int)\n",
    "        loss_pde2 = self.model2.loss_pde(x2_int)\n",
    "        loss_pde = loss_pde1 + loss_pde2\n",
    "        return loss_pde\n",
    "    #loss_ic(self, x_ic, U_ic):\n",
    "    def loss_ic(self, x_ic, U_ic):\n",
    "        x_predict = torch.zeros_like(x_ic[:, 1:2])#初始间断位置\n",
    "        \n",
    "        index1 = [i for i in range(len(x_ic)) if x_ic[i, 1:2] < x_predict[i]]\n",
    "        index2 = [i for i in range(len(x_ic)) if x_ic[i, 1:2] >= x_predict[i]]        \n",
    "        x1_ic = x_ic[index1]\n",
    "        x2_ic = x_ic[index2]\n",
    "        U1_ic = U_ic[index1]\n",
    "        U2_ic = U_ic[index2]\n",
    "        loss_ic1 = self.model1.loss_ic(x1_ic, U1_ic)\n",
    "        loss_ic2 = self.model2.loss_ic(x2_ic, U2_ic)\n",
    "        loss_ic = loss_ic1 + loss_ic2\n",
    "        return loss_ic\n",
    "    def loss_bc(self, x_leftbc, x_rightbc, U_leftbc, U_rightbc):\n",
    "        loss_bc1 = self.model1.loss_bc(x_leftbc, U_leftbc)\n",
    "        loss_bc2 = self.model2.loss_bc(x_rightbc, U_rightbc)\n",
    "        loss_bc = loss_bc1 + loss_bc2\n",
    "        return loss_bc        \n",
    "    def loss_xic(self, t_ic, x_ic):\n",
    "        x_predict = self.model3.forword(t_ic)[:,0:1]\n",
    "        loss_ic = ((x_predict - x_ic[:, 0:1]) ** 2).mean()\n",
    "        return loss_ic\n",
    "    #loss_bc(self, x_bc, U_bc):\n",
    "    #R-H，D=u关系\n",
    "    def loss_rh(self, t, x_start, x_end):\n",
    "        x_predict = self.model3.forword(t)[:,0:1]\n",
    "        dx = gradients(x_predict, t)[0]\n",
    "        dx_t = dx[:, 0:1]\n",
    "        index = torch.where((x_predict<=x_end)&(x_predict>=x_start))\n",
    "        t = t[index]\n",
    "        x_predict = x_predict[index]\n",
    "        dx_t = dx_t[index]\n",
    "        #\n",
    "        x_int = to_numpy(torch.cat([t.reshape(-1,1), x_predict.reshape(-1,1)], dim=1))\n",
    "        x_int = torch.tensor(x_int, dtype=dtype, requires_grad=True, device=device)\n",
    "        rho1, u1, p1 = self.model1.forword(x_int)\n",
    "        rho2, u2, p2 = self.model2.forword(x_int)\n",
    "        \n",
    "        loss_rh = ((u1 - u2) ** 2).mean()+ ((p1 - p2) ** 2).mean()+((dx_t - u1) ** 2).mean()+((dx_t - u2) ** 2).mean()\n",
    "        return loss_rh\n",
    "       "
   ],
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-09-24T08:23:46.105022Z",
     "iopub.execute_input": "2024-09-24T08:23:46.105460Z",
     "iopub.status.idle": "2024-09-24T08:23:46.123419Z",
     "shell.execute_reply.started": "2024-09-24T08:23:46.105418Z",
     "shell.execute_reply": "2024-09-24T08:23:46.122561Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-12-10T01:45:43.148807Z",
     "start_time": "2024-12-10T01:45:43.134196Z"
    }
   },
   "id": "6be409ae6f6ebf92",
   "execution_count": 7,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "#网格构建\n",
    "num_x = 201\n",
    "num_t = 51\n",
    "Tstart = 0.0\n",
    "Tend = 0.5\n",
    "Xstart = -1.0\n",
    "Xend = 1.0\n",
    "x = np.around(np.linspace(Xstart, Xend, num_x), decimals=3)\n",
    "t = np.around(np.linspace(Tstart, Tend, num_t), decimals=3)\n",
    "t_grid, x_grid = np.meshgrid(t, x)\n",
    "T = t_grid.flatten()[:, None]\n",
    "X = x_grid.flatten()[:, None]\n",
    "X_star = np.hstack((T.flatten()[:, None], X.flatten()[:, None]))\n",
    "# 时空边界\n",
    "lb = X_star.min(0)  # 时空下界，\n",
    "ub = X_star.max(0)  # 时空上界\n",
    "print('lb is :',lb, ', ub is :', ub)"
   ],
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-09-24T08:23:46.413829Z",
     "iopub.execute_input": "2024-09-24T08:23:46.414437Z",
     "iopub.status.idle": "2024-09-24T08:23:46.432210Z",
     "shell.execute_reply.started": "2024-09-24T08:23:46.414400Z",
     "shell.execute_reply": "2024-09-24T08:23:46.431284Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-12-10T01:45:43.164655Z",
     "start_time": "2024-12-10T01:45:43.149810Z"
    }
   },
   "id": "84aca8ea777720d8",
   "execution_count": 8,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "#初值\n",
    "#loss_ic(self, x_ic, U_ic):\n",
    "num_i_train = 201\n",
    "id_ic = np.random.choice(num_x, num_i_train, replace=False)\n",
    "irhoL, iuL, ipL = 10.0, 1.0, 0.2\n",
    "irhoR, iuR, ipR = 0.1, 1.0, 0.2\n",
    "rho_u_p_LR = [irhoL, iuL, ipL, irhoR, iuR, ipR]\n",
    "x_ic = x_grid[id_ic, 0][:, None]\n",
    "t_ic = t_grid[id_ic, 0][:, None]\n",
    "x_ic_train = np.hstack((t_ic, x_ic))\n",
    "U_ic = IC(x_ic_train, irhoL, iuL, ipL, irhoR, iuR, ipR)\n",
    "#tensor\n",
    "x_ic_train = torch.tensor(x_ic_train, dtype=dtype, requires_grad=True, device=device)\n",
    "U_ic = torch.tensor(U_ic, dtype=dtype, device=device)\n",
    "#loss_xic(self, t_ic, x_ic):\n",
    "t_ic_train = t_grid[id_ic, 0][:, None]\n",
    "x_ic = t_ic_train*1\n",
    "#tensor\n",
    "t_ic_train = torch.tensor(t_ic_train, dtype=dtype, device=device, requires_grad=True)\n",
    "x_ic = torch.tensor(x_ic, dtype=dtype, device=device)"
   ],
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-09-24T08:23:46.980811Z",
     "iopub.execute_input": "2024-09-24T08:23:46.981641Z",
     "iopub.status.idle": "2024-09-24T08:23:47.183851Z",
     "shell.execute_reply.started": "2024-09-24T08:23:46.981604Z",
     "shell.execute_reply": "2024-09-24T08:23:47.183052Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-12-10T01:46:08.815917Z",
     "start_time": "2024-12-10T01:46:08.715653Z"
    }
   },
   "id": "bb1f3322bfae61ab",
   "execution_count": 9,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-12-10T01:52:10.979349Z",
     "start_time": "2024-12-10T01:52:10.971328Z"
    }
   },
   "cell_type": "code",
   "source": [
    "#bc\n",
    "# 边界设置,all time and x=lb,x=ub\n",
    "x_lbc = x_grid[0, :][:, None]\n",
    "x_rbc = x_grid[-1, :][:, None]\n",
    "t_bc = t_grid[0, :][:, None]\n",
    "x_lbc_train = np.hstack([t_bc, x_lbc])\n",
    "x_rbc_train = np.hstack([t_bc, x_rbc])\n",
    "U_lbc = IC(x_lbc_train, irhoL, iuL, ipL, irhoR, iuR, ipR)\n",
    "U_rbc = IC(x_rbc_train, irhoL, iuL, ipL, irhoR, iuR, ipR)\n",
    "#\n",
    "x_lbc_train = torch.tensor(x_lbc_train, dtype=dtype, requires_grad=True, device=device)\n",
    "U_lbc = torch.tensor(U_lbc, dtype=dtype, device=device)\n",
    "x_rbc_train = torch.tensor(x_rbc_train, dtype=dtype, requires_grad=True, device=device)\n",
    "U_rbc = torch.tensor(U_rbc, dtype=dtype, device=device)"
   ],
   "id": "da67df7f1bb51382",
   "execution_count": 10,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": "#残差\n#loss_pde(self, x_int):\nnum_f_train = 5000\nid_f = np.random.choice(num_x * num_t, num_f_train, replace=False)\nx_int = X[:, 0][id_f, None]\nt_int = T[:, 0][id_f, None]\nx_int_train = np.hstack((t_int, x_int))\n#tensor\nx_int_train = torch.tensor(x_int_train, dtype=dtype, requires_grad=True, device=device)",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-09-24T08:23:47.599759Z",
     "iopub.execute_input": "2024-09-24T08:23:47.600812Z",
     "iopub.status.idle": "2024-09-24T08:23:47.609794Z",
     "shell.execute_reply.started": "2024-09-24T08:23:47.600765Z",
     "shell.execute_reply": "2024-09-24T08:23:47.608821Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-12-10T01:52:12.704645Z",
     "start_time": "2024-12-10T01:52:12.692942Z"
    }
   },
   "id": "b5ad1a076b1101fa",
   "execution_count": 11,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "#rh关系\n",
    "#loss_rh(self, t, x_start, x_end)\n",
    "t_all = t_grid[0, :][:, None]\n",
    "x_start = Xstart\n",
    "x_end = Xend\n",
    "x_exact = t_all *1.0\n",
    "#tensor\n",
    "t_all = torch.tensor(t_all, dtype=dtype, requires_grad=True, device=device)\n",
    "x_start = torch.tensor(x_start, dtype=dtype, device=device)\n",
    "x_end = torch.tensor(x_end, dtype=dtype, device=device)"
   ],
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-09-24T08:23:48.313405Z",
     "iopub.execute_input": "2024-09-24T08:23:48.314162Z",
     "iopub.status.idle": "2024-09-24T08:23:48.321599Z",
     "shell.execute_reply.started": "2024-09-24T08:23:48.314120Z",
     "shell.execute_reply": "2024-09-24T08:23:48.320665Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-12-10T01:53:38.675546Z",
     "start_time": "2024-12-10T01:53:38.666523Z"
    }
   },
   "id": "b889c0c6fa2ecdf9",
   "execution_count": 14,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "#实例化神经网络\n",
    "layers1 = Dnnlayers(2, 3, 3, 10)\n",
    "layers2 = Dnnlayers(2, 3, 3, 10)\n",
    "layers3 = Dnnlayers(1, 1, 2, 5)\n",
    "# Pinns(layers)\n",
    "model = conbine(layers1, layers2, layers3)\n",
    "loss_total_history = []\n",
    "loss_pde_history = []\n",
    "loss_ibc_history = []\n",
    "loss_rh_history = []\n",
    "lr = 0.001\n",
    "optimizer = torch.optim.Adam(list(model.model1.dnn.parameters())+list(model.model2.dnn.parameters())+list(model.model3.parameters()), lr=lr)"
   ],
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-09-24T08:23:49.532299Z",
     "iopub.execute_input": "2024-09-24T08:23:49.533197Z",
     "iopub.status.idle": "2024-09-24T08:23:50.456244Z",
     "shell.execute_reply.started": "2024-09-24T08:23:49.533152Z",
     "shell.execute_reply": "2024-09-24T08:23:50.455427Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-12-10T01:53:40.675325Z",
     "start_time": "2024-12-10T01:53:40.666301Z"
    }
   },
   "id": "70089af1370df782",
   "execution_count": 15,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "#training\n",
    "def train(epoch):\n",
    "    model.model1.dnn.train()\n",
    "    model.model2.dnn.train()\n",
    "    model.model3.train()\n",
    "    def closure():\n",
    "        optimizer.zero_grad()\n",
    "\n",
    "        loss_pde = model.loss_pde(x_int_train)\n",
    "        loss_ic = model.loss_ic(x_ic_train, U_ic)\n",
    "        loss_bc = model.loss_bc(x_lbc_train, x_rbc_train, U_lbc, U_rbc)\n",
    "        loss_xic = model.loss_xic(t_ic_train, x_ic)\n",
    "        loss_rh = model.loss_rh(t_all, x_start, x_end)\n",
    "                   \n",
    "        loss = loss_pde + 10*(loss_ic+loss_bc + loss_rh + loss_xic) \n",
    "        if epoch%1==0:\n",
    "            print(f'epoch:{epoch}, loss:{loss:.8f},loss_pde:{loss_pde:.8f}, loss_IBC:{loss_ic+loss_bc:.8f}, loss_xic:{loss_xic:.8f},loss_rh:{loss_rh:.8f}')\n",
    "        if epoch<300:\n",
    "            if epoch%10==0:\n",
    "                x_predict = to_numpy(model.model3.forword(t_all)[:,0:1])\n",
    "                plt.figure()\n",
    "                plt.plot(x_exact, to_numpy(t_all), c = 'r', label='x_exact')\n",
    "                plt.scatter(x_predict, to_numpy(t_all), cmap='viridis', alpha=0.8, label='x_predict')\n",
    "                plt.xlabel('x')\n",
    "                plt.ylabel('t')\n",
    "                plt.title(f'epoch is {epoch}')\n",
    "                plt.legend()\n",
    "                plt.savefig(f'the {epoch} fitting of x(t).png')\n",
    "                plt.show()\n",
    "        if epoch%1000==0:\n",
    "            x_predict = to_numpy(model.model3.forword(t_all)[:,0:1])\n",
    "            plt.figure()\n",
    "            plt.plot(x_exact, to_numpy(t_all), c = 'r', label='x_exact')\n",
    "            plt.scatter(x_predict, to_numpy(t_all), cmap='viridis', alpha=0.8, label='x_predict')\n",
    "            plt.xlabel('x')\n",
    "            plt.ylabel('t')\n",
    "            plt.title(f'epoch is {epoch}')\n",
    "            plt.legend()\n",
    "            plt.savefig(f'the {epoch} fitting of x(t).png')\n",
    "            plt.show()\n",
    "\n",
    "\n",
    "        loss_total_history.append(to_numpy(loss))\n",
    "        loss_pde_history.append(to_numpy(loss_pde))\n",
    "        loss_ibc_history.append(to_numpy(loss_ic+loss_xic+loss_bc))\n",
    "        loss_rh_history.append(to_numpy(loss_rh))\n",
    "        \n",
    "        loss.backward()\n",
    "        return loss\n",
    "    loss = optimizer.step(closure)"
   ],
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-09-24T08:23:50.457846Z",
     "iopub.execute_input": "2024-09-24T08:23:50.458658Z",
     "iopub.status.idle": "2024-09-24T08:23:50.466498Z",
     "shell.execute_reply.started": "2024-09-24T08:23:50.458610Z",
     "shell.execute_reply": "2024-09-24T08:23:50.465527Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-12-10T01:55:51.195182Z",
     "start_time": "2024-12-10T01:55:51.172376Z"
    }
   },
   "id": "f674917808bcff0d",
   "execution_count": 19,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-12-10T01:55:52.116325Z",
     "start_time": "2024-12-10T01:55:51.945664Z"
    }
   },
   "cell_type": "code",
   "source": [
    "x_predict = to_numpy(model.model3.forword(t_all)[:,0:1])\n",
    "plt.figure()\n",
    "plt.plot(x_exact, to_numpy(t_all), c = 'r', label='x_exact')\n",
    "plt.scatter(x_predict, to_numpy(t_all), cmap='viridis', alpha=0.8, label='x_predict')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title(f'epoch is {0}')\n",
    "plt.legend()\n",
    "plt.savefig(f'the {0} fitting of x(t).png')\n",
    "plt.show()\n"
   ],
   "id": "43672f8fec9fb9d2",
   "execution_count": 20,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": "epochs1 = 10000\n#for training\ntic = time.time()\nfor epoch in range(1, epochs1+1):\n    train(epoch)\ntoc = time.time()",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-09-24T08:23:52.150100Z",
     "iopub.execute_input": "2024-09-24T08:23:52.151126Z",
     "iopub.status.idle": "2024-09-24T09:31:41.292026Z",
     "shell.execute_reply.started": "2024-09-24T08:23:52.151082Z",
     "shell.execute_reply": "2024-09-24T09:31:41.291134Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-12-10T02:09:10.168537Z",
     "start_time": "2024-12-10T01:55:53.772402Z"
    }
   },
   "id": "59e6753e762b2257",
   "execution_count": 21,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "optimizer = torch.optim.LBFGS(list(model.model1.dnn.parameters())+list(model.model2.dnn.parameters())+list(model.model3.parameters()),lr=1.0,max_iter=10)\n",
    "epochs2 = 1000\n",
    "tic = time.time()\n",
    "for epoch in range(15001, 15000+epochs2+1):\n",
    "    train(epoch)\n",
    "    if (epoch-15000)%100==0:\n",
    "        x_predict = to_numpy(model.model3.forword(t_all)[:,0:1])\n",
    "        plt.figure()\n",
    "        plt.plot(x_exact, to_numpy(t_all), c = 'r', label='x_exact')\n",
    "        plt.scatter(x_predict, to_numpy(t_all), cmap='viridis', alpha=0.8, label='x_predict')\n",
    "        plt.xlabel('x')\n",
    "        plt.ylabel('t')\n",
    "        plt.title(f'epoch is {epoch}')\n",
    "        plt.legend()\n",
    "        plt.savefig(f'the {epoch} fitting of x(t).png')\n",
    "        plt.show()\n",
    "toc = time.time()"
   ],
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-09-24T09:31:41.293728Z",
     "iopub.execute_input": "2024-09-24T09:31:41.294028Z",
     "iopub.status.idle": "2024-09-24T09:38:28.107433Z",
     "shell.execute_reply.started": "2024-09-24T09:31:41.293996Z",
     "shell.execute_reply": "2024-09-24T09:38:28.106494Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-12-10T02:22:20.116598Z",
     "start_time": "2024-12-10T02:20:47.974402Z"
    }
   },
   "id": "cdd6739e70b2ed04",
   "execution_count": 23,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-12-10T02:22:20.585594Z",
     "start_time": "2024-12-10T02:22:20.118603Z"
    }
   },
   "cell_type": "code",
   "source": [
    "eporchx = np.arange(1, len(loss_total_history)+1)\n",
    "plt.plot(eporchx, loss_total_history, label='loss_total_history')\n",
    "plt.plot(eporchx, loss_pde_history, label='loss_pde_history')\n",
    "plt.plot(eporchx, loss_ibc_history, label='loss_ibc_history')\n",
    "plt.plot(eporchx, loss_rh_history, label='loss_rh_history')\n",
    "plt.yscale('log')\n",
    "plt.title('loss with eproch')\n",
    "plt.legend()\n",
    "plt.show()"
   ],
   "id": "5d53fe4a435204b4",
   "execution_count": 24,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-12-10T02:22:21.601398Z",
     "start_time": "2024-12-10T02:22:20.586912Z"
    }
   },
   "cell_type": "code",
   "source": [
    "x_int_test = np.hstack([T[:, 0][:, None], X[:, 0][:, None]])\n",
    "\n",
    "#exact\n",
    "rho_exact = np.zeros_like(x_int_test[:, 0:1])\n",
    "u_exact = np.zeros_like(x_int_test[:, 0:1])\n",
    "p_exact = np.zeros_like(x_int_test[:, 0:1])\n",
    "u_exact[:, 0:1] = 1.0\n",
    "p_exact[:, 0:1] = 0.2\n",
    "x_exact = to_numpy(x_int_test[:, 0:1])*1.0\n",
    "\n",
    "index = np.where(x_int_test[:, 1:2]<x_exact[:,0:1])\n",
    "print('index:', len(index))\n",
    "index1 = [i for i in range(len(x_int_test)) if x_int_test[i, 1:2] < x_exact[i]]\n",
    "index2 = [i for i in range(len(x_int_test)) if x_int_test[i, 1:2] >= x_exact[i]]\n",
    "rho_exact[index1] = 10.0\n",
    "rho_exact[index2] = 0.1\n",
    "x_int_test1 = x_int_test[index1]\n",
    "index = np.where(x_int_test1[:, 1]==x_int_test1[:, 0])\n",
    "print('index1:', len(index1))\n",
    "print('index:', index)\n",
    "t_grid_test = to_numpy(x_int_test[:, 0:1]).reshape(num_x, num_t)\n",
    "x_grid_test = to_numpy(x_int_test[:, 1:2]).reshape(num_x, num_t)\n",
    "rho_exact = rho_exact.reshape(num_x, num_t)\n",
    "u_exact = u_exact.reshape(num_x, num_t)\n",
    "p_exact = p_exact.reshape(num_x, num_t)\n",
    "scipy.io.savemat(f'1Dexact_Case1.mat', {'x': x, 't': t, 'rho': rho_exact.T, 'u': u_exact.T, 'p': p_exact.T})\n",
    "\n",
    "plt.figure()\n",
    "plt.scatter(x_int_test[:, 1:2], x_int_test[:, 0:1], cmap='viridis', alpha=1, s=0.005)\n",
    "plt.scatter(x_exact, x_int_test[:, 0:1], cmap='viridis', alpha=1.0, s=0.005)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.show()\n",
    "plt.figure()\n",
    "plt.scatter(x_int_test1[:, 1:2], x_int_test1[:, 0:1], cmap='viridis', alpha=1.0, s=0.005)\n",
    "plt.scatter(x_exact, x_int_test[:, 0:1], cmap='viridis', alpha=1.0, s=0.005)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.show()\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, rho_exact, levels=np.linspace(rho_exact.min(), rho_exact.max(), 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('rho_exact')\n",
    "plt.show()"
   ],
   "id": "9cad8285ce5d7ca9",
   "execution_count": 25,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "t_test = to_numpy(x_int_test[:, 0:1])\n",
    "t_test = torch.tensor(t_test, dtype=dtype, requires_grad=True, device=device)\n",
    "x_test = x_int_test[:, 1:2]\n",
    "x_predict = np.around(to_numpy(model.model3.forword(t_test)[:,0:1]), decimals=3)\n",
    "plt.figure()\n",
    "plt.scatter(x_predict, to_numpy(t_test), cmap='viridis', alpha=0.8, label='x_predict', s=0.005)\n",
    "#plt.plot(x_exact, to_numpy(t_test), 'b', label='x_exact')\n",
    "plt.scatter(x_exact, to_numpy(t_test), label='x_exact', s=0.005)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "error = abs(x_predict - x_exact)\n",
    "plt.figure()\n",
    "plt.scatter(to_numpy(t_test), error, cmap='viridis', alpha=0.8, label='x_error')\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('error')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "index1 = [i for i in range(len(x_int_test)) if x_test[i] < x_predict[i]]\n",
    "index2 = [i for i in range(len(x_int_test)) if x_test[i] >= x_predict[i]]\n",
    "\n",
    "x1_int = x_int_test[index1]\n",
    "'''x_predict1 = x_predict[index1]\n",
    "print(len(index1))\n",
    "print('x_predict1[5100]',x_predict1[5100])\n",
    "mask = np.all(np.any(x1_int[:, None] != x_int_test1, axis=2), axis=1)\n",
    "print([i for i in range(len(mask)) if mask[i]])\n",
    "cha = x1_int[mask]\n",
    "print('index1:', len(index1))\n",
    "print(cha)\n",
    "plt.figure()\n",
    "plt.scatter(x_int_test1[:, 1:2], x_int_test1[:, 0:1], cmap='viridis', alpha=1.0, s=0.005)\n",
    "plt.scatter(x1_int[:, 1:2], x1_int[:, 0:1], cmap='viridis', alpha=1.0, s=0.005)\n",
    "#plt.scatter(x_exact, x_int_test[:, 0:1], cmap='viridis', alpha=1.0, s=0.005)\n",
    "#plt.scatter(x_predict, x_int_test[:, 0:1], cmap='viridis', alpha=1.0, s=0.005)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.show()'''\n",
    "\n",
    "x2_int = x_int_test[index2]\n",
    "x1_int = torch.tensor(x1_int, dtype=dtype, requires_grad=True, device=device)\n",
    "x2_int = torch.tensor(x2_int, dtype=dtype, requires_grad=True, device=device)\n",
    "rho_predict1, u_predict1, p_predict1 = model.model1.forword(x1_int)\n",
    "rho_predict2, u_predict2, p_predict2 = model.model2.forword(x2_int)\n",
    "\n",
    "rho_predict = torch.vstack([rho_predict1, rho_predict2])\n",
    "rho_predict[index1] = rho_predict1\n",
    "rho_predict[index2] = rho_predict2\n",
    "u_predict = torch.vstack([u_predict1, u_predict2])\n",
    "u_predict[index1] = u_predict1\n",
    "u_predict[index2] = u_predict2\n",
    "p_predict = torch.vstack([p_predict1, p_predict2])\n",
    "p_predict[index1] = p_predict1\n",
    "p_predict[index2] = p_predict2\n",
    "\n",
    "\n",
    "rho_predict = to_numpy(rho_predict).reshape(num_x, num_t)\n",
    "u_predict = to_numpy(u_predict).reshape(num_x, num_t)\n",
    "p_predict = to_numpy(p_predict).reshape(num_x, num_t)\n",
    "\n",
    "\n",
    "\n",
    "scipy.io.savemat(f'1Dpredict_Case1.mat', {'x': x, 't': t, 'rho': rho_predict.T, 'u': u_predict.T, 'p': p_predict.T})\n"
   ],
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-12-10T02:22:22.256803Z",
     "start_time": "2024-12-10T02:22:21.603405Z"
    }
   },
   "id": "f56e33350875c610",
   "execution_count": 26,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-12-10T02:22:26.723128Z",
     "start_time": "2024-12-10T02:22:22.257806Z"
    }
   },
   "cell_type": "code",
   "source": [
    "\n",
    "error_rho = abs(rho_predict - rho_exact)\n",
    "error_u = abs(u_predict - u_exact)\n",
    "error_p = abs(p_predict - p_exact)\n",
    "#\n",
    "\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, rho_predict, levels=np.linspace(rho_predict.min(), rho_predict.max(), 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('rho_predict')\n",
    "plt.show()\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, u_predict, levels=np.linspace(u_predict.min(), u_predict.max(), 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('u_predict')\n",
    "plt.show()\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, p_predict, levels=np.linspace(p_predict.min(), p_predict.max(), 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('p_predict')\n",
    "plt.show()\n",
    "#\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, error_rho, levels=np.linspace(error_rho.min(), error_rho.max()+0.05, 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('error_rho')\n",
    "plt.show()\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, error_u, levels=np.linspace(error_u.min(), error_u.max(), 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('error_u')\n",
    "plt.show()\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, error_p, levels=np.linspace(error_p.min(), error_p.max(), 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('error_p')\n",
    "plt.show()\n",
    "\n",
    "print('L2 for error_rho is :', np.sqrt(((error_rho)**2).mean()))\n",
    "print('L2 for error_u is :', np.sqrt(((error_u)**2).mean()))\n",
    "print('L2 for error_p is :', np.sqrt(((error_p)**2).mean()))\n"
   ],
   "id": "950ca570707fd5ed",
   "execution_count": 27,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-12-10T06:13:57.795706Z",
     "start_time": "2024-12-10T06:13:49.011716Z"
    }
   },
   "cell_type": "code",
   "source": [
    "data = scipy.io.loadmat(r'1Dexact_Case1.mat')  # Import Solution data\n",
    "x_test = data['x'].flatten()[:, None]  # Partitioned spatial coordinates\n",
    "t_test = data['t'].flatten()[:, None]  # Partitioned spatial coordinates\n",
    "rho_exact = np.real(data['rho'])  # rho(t,x,y)\n",
    "p_exact = np.real(data['p'])\n",
    "u_exact = np.real(data['u'])  # Exact total\n",
    "data = scipy.io.loadmat(r'1Dpredict_Case1.mat')  # Import Solution data\n",
    "rho_predict = np.real(data['rho'])\n",
    "p_predict = np.real(data['p'])\n",
    "u_predict = np.real(data['u'])\n",
    "\n",
    "t_grid_test, x_grid_test = np.meshgrid(t_test, x_test)\n",
    "x_int_test = np.hstack((t_grid_test.flatten()[:, None], x_grid_test.flatten()[:, None]))\n",
    "x_int_test = torch.tensor(x_int_test, dtype=dtype, requires_grad=True, device=device)\n",
    "\n",
    "\n",
    "\n",
    "error_rho = abs(rho_predict - rho_exact)\n",
    "error_u = abs(u_predict - u_exact)\n",
    "error_p = abs(p_predict - p_exact)\n",
    "#\n",
    "\n",
    "plt.figure()\n",
    "c = plt.contourf(x_grid_test, t_grid_test, rho_exact.T, levels=np.linspace(rho_exact.min(), rho_exact.max(), 1000), cmap='jet')\n",
    "# 添加颜色条，并设置刻度范围\n",
    "ticks = [0.1, 2, 4, 6, 8, 10]  # 自定义刻度\n",
    "cbar = plt.colorbar(c, ticks=ticks)\n",
    "cbar.set_ticklabels(ticks)  # 设置刻度标签\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('rho_exact')\n",
    "plt.savefig('1Drho_exact_Case1.png')\n",
    "plt.show()\n",
    "\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, u_exact.T, levels=np.linspace(u_exact.min(), u_exact.max()+0.05, 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('u_exact')\n",
    "plt.savefig('1Du_exact_Case1.png')\n",
    "plt.show()\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, p_exact.T, levels=np.linspace(p_exact.min(), p_exact.max()+0.05, 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('p_exact')\n",
    "plt.savefig('1Dp_exact_Case1.png')\n",
    "plt.show()\n",
    "\n",
    "\n",
    "plt.figure()\n",
    "c = plt.contourf(x_grid_test, t_grid_test, rho_predict.T, levels=np.linspace(rho_predict.min(), rho_predict.max(), 1000), cmap='jet')\n",
    "# 添加颜色条，并设置刻度范围\n",
    "ticks = [0.1, 2, 4, 6, 8, 10]  # 自定义刻度\n",
    "cbar = plt.colorbar(c, ticks=ticks)\n",
    "cbar.set_ticklabels(ticks)  # 设置刻度标签\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('rho_predict')\n",
    "plt.savefig('1Drho_predict_Case1.png')\n",
    "\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, u_predict.T, levels=np.linspace(u_predict.min(), u_predict.max(), 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('u_predict')\n",
    "plt.savefig('1Du_predict_Case1.png')\n",
    "plt.show()\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, p_predict.T, levels=np.linspace(p_predict.min(), p_predict.max(), 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('p_predict')\n",
    "plt.savefig('1Dp_predict_Case1.png')\n",
    "plt.show()\n",
    "#\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, error_rho.T, levels=np.linspace(error_u.min(), error_u.max(), 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('error_rho')\n",
    "plt.savefig('1Derror_rho_Case1.png')\n",
    "plt.show()\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, error_u.T, levels=np.linspace(error_u.min(), error_u.max(), 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('error_u')\n",
    "plt.savefig('1Derror_u_Case1.png')\n",
    "plt.show()\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, error_p.T, levels=np.linspace(error_p.min(), error_p.max(), 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('error_p')\n",
    "plt.savefig('1Derror_p_Case1.png')\n",
    "plt.show()\n",
    "\n",
    "print('L2 for error_rho is :', np.sqrt(((error_rho)**2).mean()))\n",
    "print('L2 for error_u is :', np.sqrt(((error_u)**2).mean()))\n",
    "print('L2 for error_p is :', np.sqrt(((error_p)**2).mean()))\n",
    "print(error_rho.max())\n",
    "'''\n",
    "L2 for error_rho is : 0.00021239900986887226\n",
    "L2 for error_u is : 0.0002779203887260891\n",
    "L2 for error_p is : 8.534069059495541e-05\n",
    "'''"
   ],
   "id": "cfe3b3b1a7475a2c",
   "execution_count": 36,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-12-10T02:22:36.039190Z",
     "start_time": "2024-12-10T02:22:35.835161Z"
    }
   },
   "cell_type": "code",
   "source": [
    "plt.figure()\n",
    "plt.plot(x_test, rho_exact[-1,:], 'b', label='end_rho_exact')\n",
    "plt.scatter(x_test, rho_predict[-1,:], color='r', s=1, label='end_rho_predict')\n",
    "plt.legend()\n",
    "#plt.savefig('1Drho_end_Case1.png')\n",
    "plt.show()\n",
    "plt.figure()\n",
    "plt.plot(x_test[10:150], rho_exact[-1,10:150], 'b', label='end_rho_exact')\n",
    "plt.scatter(x_test[10:150], rho_predict[-1,10:150], color='r', s=1, label='end_rho_predict')\n",
    "plt.legend()\n",
    "#plt.savefig('1Drho_end_Case1.png')\n",
    "plt.show()"
   ],
   "id": "83fe94d235434707",
   "execution_count": 29,
   "outputs": []
  },
  {
   "metadata": {},
   "cell_type": "code",
   "execution_count": null,
   "source": "",
   "id": "b25b7fedbf47fc7d",
   "outputs": []
  }
 ]
}
